#include <cstdio>
#include <cstring>
#include <iostream>
#include "Riostream.h"
#include "TROOT.h"
#include "TGraph.h"

#include "SRIMloader.hh"

ifstream SRIMfile;

void SRIM::Creategraphs(char* dir)
{
	this->dir = dir;

//	SRIM::SRIMreader("he6pe"); // First time reader doesn't work!!
	// ------------- Deadlayer ------------

	// Protons
	SRIM::SRIMreader("h1si");
	TGraph *RangeDLpp = new TGraph(131,Eion,ProjRange);
	RDL[1] = RangeDLpp;
	TGraph *EnergyDLpp = new TGraph(131,ProjRange,Eion);
	EDL[1] = EnergyDLpp;

	// Deuterium
	SRIM::SRIMreader("h2si");
	TGraph *RangeDLdp = new TGraph(131,Eion,ProjRange);
	RDL[2] = RangeDLdp;
	TGraph *EnergyDLdp = new TGraph(131,ProjRange,Eion);
	EDL[2] = EnergyDLdp;

	// Tritium
	SRIM::SRIMreader("h3si");
	TGraph *RangeDLtp = new TGraph(131,Eion,ProjRange);
	RDL[3] = RangeDLtp;
	TGraph *EnergyDLtp = new TGraph(131,ProjRange,Eion);
	EDL[3] = EnergyDLtp;

	// Alpha
	SRIM::SRIMreader("he4si");
	TGraph *RangeDLap = new TGraph(131,Eion,ProjRange);
	RDL[4] = RangeDLap;
	TGraph *EnergyDLap = new TGraph(131,ProjRange,Eion);
	EDL[4] = EnergyDLap;

	// Beryllium 12
	SRIM::SRIMreader("be12si");
	TGraph *RangeDL12bep = new TGraph(131,Eion,ProjRange);
	RDL[12] = RangeDL12bep;
	TGraph *EnergyDL12bep = new TGraph(131,ProjRange,Eion);
	EDL[12] = EnergyDL12bep;

	// Beryllium 11
	SRIM::SRIMreader("be11si");
	TGraph *RangeDL11bep = new TGraph(131,Eion,ProjRange);
	RDL[11] = RangeDL11bep;
	TGraph *EnergyDL11bep = new TGraph(131,ProjRange,Eion);
	EDL[11] = EnergyDL11bep;

	// Beryllium 10
	SRIM::SRIMreader("be10si");
	TGraph *RangeDL10bep = new TGraph(131,Eion,ProjRange);
	RDL[10] = RangeDL10bep;
	TGraph *EnergyDL10bep = new TGraph(131,ProjRange,Eion);
	EDL[10] = EnergyDL10bep;
	
	// Lithium 9
	SRIM::SRIMreader("li9si");
	TGraph *RangeDL9lip = new TGraph(131,Eion,ProjRange);
	RDL[9] = RangeDL9lip;
	TGraph *EnergyDL9lip = new TGraph(131,ProjRange,Eion);
	EDL[9] = EnergyDL9lip;
	
	// Carbon 9
	SRIM::SRIMreader("c12si");
	TGraph *RangeDL12cp = new TGraph(131,Eion,ProjRange);
	RDL[13] = RangeDL12cp;
	TGraph *EnergyDL12cp = new TGraph(131,ProjRange,Eion);
	EDL[13] = EnergyDL12cp;
	
	// ------------------- Target -------------------

	// Protons
	SRIM::SRIMreader("h1pe");
	TGraph *RangeTpp = new TGraph(131,Eion,ProjRange);
	RT[1] = RangeTpp;
	TGraph *EnergyTpp = new TGraph(131,ProjRange,Eion);
	ET[1] = EnergyTpp;

	// Alpha
	SRIM::SRIMreader("he4pe");
	TGraph *RangeTap = new TGraph(131,Eion,ProjRange);
	RT[4] = RangeTap;
	TGraph *EnergyTap = new TGraph(131,ProjRange,Eion);
	ET[4] = EnergyTap;

	// Tritium
	SRIM::SRIMreader("h3pe");
	TGraph *RangeTtp = new TGraph(131,Eion,ProjRange);
	RT[3] = RangeTtp;
	TGraph *EnergyTtp = new TGraph(131,ProjRange,Eion);
	ET[3] = EnergyTtp;
	
	// Beryllium 12
	SRIM::SRIMreader("be12pe");
	TGraph *RangeT12bep = new TGraph(131,Eion,ProjRange);
	RT[12] = RangeT12bep;
	TGraph *EnergyT12bep = new TGraph(131,ProjRange,Eion);
	ET[12] = EnergyT12bep;
	
	// Deuterium
	SRIM::SRIMreader("h2pe");
	TGraph *RangeTdp = new TGraph(131,Eion,ProjRange);
	RT[2] = RangeTdp;
	TGraph *EnergyTdp = new TGraph(131,ProjRange,Eion);
	ET[2] = EnergyTdp;
	
	// Beryllium 11
	SRIM::SRIMreader("be11pe");
	TGraph *RangeT11bep = new TGraph(131,Eion,ProjRange);
	RT[11] = RangeT11bep;
	TGraph *EnergyT11bep = new TGraph(131,ProjRange,Eion);
	ET[11] = EnergyT11bep;
	
	// Beryllium 10
	SRIM::SRIMreader("be10pe");
	TGraph *RangeT10bep = new TGraph(131,Eion,ProjRange);
	RT[10] = RangeT10bep;
	TGraph *EnergyT10bep = new TGraph(131,ProjRange,Eion);
	ET[10] = EnergyT10bep;
/*	
	// Lithium 9
	SRIM::SRIMreader("li9pe");
	TGraph *RangeT9lip = new TGraph(131,Eion,ProjRange);
	RT[9] = RangeT9lip;
	TGraph *EnergyT9lip = new TGraph(131,ProjRange,Eion);
	ET[9] = EnergyT9lip;
*/	
	
	// Carbon 12
	SRIM::SRIMreader("c12pe");
	TGraph *RangeT12cp = new TGraph(131,Eion,ProjRange);
	RT[13] = RangeT12cp;
	TGraph *EnergyT12cp = new TGraph(131,ProjRange,Eion);
	ET[13] = EnergyT12cp;
	
	// Beryllium 11 on silver
	SRIM::SRIMreader("be11ag");
	TGraph *RangeT11beag = new TGraph(131,Eion,ProjRange);
	RT[20] = RangeT11beag;
	TGraph *EnergyT11beag = new TGraph(131,ProjRange,Eion);
	ET[20] = EnergyT11beag;
	
	// ------------------- Mylar -------------------

	// Protons
	SRIM::SRIMreader("h1mylar");
	TGraph *RangeMylarpp = new TGraph(131,Eion,ProjRange);
	RMylar[1] = RangeMylarpp;
	TGraph *EnergyMylarpp = new TGraph(131,ProjRange,Eion);
	EMylar[1] = EnergyMylarpp;

	// Alpha
	SRIM::SRIMreader("he4mylar");
	TGraph *RangeMylarap = new TGraph(131,Eion,ProjRange);
	RMylar[4] = RangeMylarap;
	TGraph *EnergyMylarap = new TGraph(131,ProjRange,Eion);
	EMylar[4] = EnergyMylarap;

	// Tritium
	SRIM::SRIMreader("h3mylar");
	TGraph *RangeMylartp = new TGraph(131,Eion,ProjRange);
	RMylar[3] = RangeMylartp;
	TGraph *EnergyMylartp = new TGraph(131,ProjRange,Eion);
	EMylar[3] = EnergyMylartp;
	
	// Beryllium 11
	SRIM::SRIMreader("be11mylar");
	TGraph *RangeMylar11bep = new TGraph(131,Eion,ProjRange);
	RMylar[11] = RangeMylar11bep;
	TGraph *EnergyMylar11bep = new TGraph(131,ProjRange,Eion);
	EMylar[11] = EnergyMylar11bep;
	
	// Deuterium
	SRIM::SRIMreader("h2mylar");
	TGraph *RangeMylardp = new TGraph(131,Eion,ProjRange);
	RMylar[2] = RangeMylardp;
	TGraph *EnergyMylardp = new TGraph(131,ProjRange,Eion);
	EMylar[2] = EnergyMylardp;
	
	return;
}


void SRIM::SRIMreader(char* filename)
{
	Float_t dEdxe[131], dEdxNuc[131], Long[131], Lateral[131];
	string EUnit = "keV";
	string RUnit = "A";
	string LongUnit = "A";
	string LateralUnit = "A";
	string line;
	Int_t iline = 1;
        std::string SRIMfilename = dir + '/' + filename;

	SRIMfile.open(SRIMfilename.c_str());

	// Loop to read input in file
	while (1)
	{
		if (!SRIMfile.good()) break;
		
		// Extra info, not used
		if (iline <= 25 || iline >= 157) getline (SRIMfile,line);

		// These lines are to be used
		if (iline > 25 && iline < 157)
		{	
			SRIMfile >> Eion[iline - 26] >> EUnit >> dEdxe[iline - 26] >> dEdxNuc[iline - 26] >> ProjRange[iline - 26] >> RUnit >> Long[iline - 26] >> LongUnit >> Lateral[iline - 26] >> LateralUnit;
			if (EUnit == "eV") Eion[iline - 26] = Eion[iline - 26]*1e-3; // Eion in keV
			if (EUnit == "MeV") Eion[iline - 26] = Eion[iline - 26]*1e+3; // Eion in keV
			if (RUnit == "A") ProjRange[iline - 26] = ProjRange[iline - 26]*1e-4; // Range in ym
			if (RUnit == "mm") ProjRange[iline - 26] = ProjRange[iline - 26]*1e+3; // Range in ym
		}
	iline++;
	}
	SRIMfile.close();


	return;
}
